Understanding the model used by lme4

lme4’s model

Simplified from Bates et al. (2015).

Linear model conditional on random effects \(b\):
\[ (y \,|\,b) \mathrel{\sim}\mathcal{N} \left( X \beta + Z b, \sigma^2 I_n \right) \]
\(X\) is an \(n \times p\) design matrix.

Distribution of \(b\):
\[ b \mathrel{\sim}\mathcal{N} (0, \Sigma) \]
\(Z\) is an \(n \times q\) design matrix, \(\Sigma\) is positive definite.

More precisely, there can be several random effects terms \(i\), \(Z_i b_i\), and for each \(Z_i\) contains a \(q_i\) regressors within \(r_i\) groups. The coefficients for the regressors within each group have a \(q_i × q_i\)
covariance \(\Sigma_i\) across groups. I.e.
\[ b_i \mathrel{\sim}\mathcal{N} (0, I_{r_i} ⊗ \Sigma_i) \]

The example data

Load sleepstudy and restrict to 3 days and 3 subjects:

library(lme4, quietly=TRUE)
data <- sleepstudy[(sleepstudy$Days <= 2) & (sleepstudy$Subject %in% c(308, 309, 310)), ]
data

	Reaction	Days	Subject
1	249.5600	0	308
2	258.7047	1	308
3	250.8006	2	308
11	222.7339	0	309
12	205.2658	1	309
13	202.9778	2	309
21	199.0539	0	310
22	194.3322	1	310
23	234.3200	2	310

Model, variables are (Intercept) and Days, both fixed and random:

model <- lmer(Reaction ~ Days + (Days | Subject), data)
summary(model)

Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (Days | Subject)
   Data: data

REML criterion at convergence: 65.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3294 -0.3361 -0.0463  0.5375  0.9859 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 Subject  (Intercept) 808.5    28.43         
          Days        122.8    11.08    -0.56
 Residual             140.0    11.83         
Number of obs: 9, groups:  Subject, 3

Fixed effects:
            Estimate Std. Error t value
(Intercept)  221.403     17.561  12.607
Days           2.792      8.016   0.348

Correlation of Fixed Effects:
     (Intr)
Days -0.585

We now try to extract the parts of the model in the equation above (constants and estimates).

fixed effects

model matrix

getME(model, "X")

   (Intercept) Days
1            1    0
2            1    1
3            1    2
11           1    0
12           1    1
13           1    2
21           1    0
22           1    1
23           1    2
attr(,"assign")
[1] 0 1
attr(,"msgScaleX")
character(0)

c(getME(model, "n"), getME(model, "p"))

[1] 9 2

parameter estimates

# getME(model, "beta")
getME(model, "fixef")

(Intercept)        Days 
 221.402556    2.791767 

random effects

model matrix

getME(model, "Z")

9 x 6 sparse Matrix of class "dgCMatrix"
   308 308 309 309 310 310
1    1   .   .   .   .   .
2    1   1   .   .   .   .
3    1   2   .   .   .   .
11   .   .   1   .   .   .
12   .   .   1   1   .   .
13   .   .   1   2   .   .
21   .   .   .   .   1   .
22   .   .   .   .   1   1
23   .   .   .   .   1   2

c(getME(model, "n"), getME(model, "q"))

[1] 9 6

residual standard deviation

# getME(model, "sigma")
sigma(model)

[1] 11.83223

getME(model, "theta")

     Subject.(Intercept) Subject.Days.(Intercept)             Subject.Days 
               2.4031809               -0.5266796                0.7742044